Entanglement of remote quantum systems by environmental modes 
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We investigate the generation of quantum mechanical entanglement of two remote osciUators that 
are locaUy coupled to a common bosonic bath. Starting with a Lagrangian formulation of a suitable 
model, we derive two coupled Quantum Langevin Equations that exactly describe the time evolution 
of the two local oscillators in presence of the coupling to the bosonic bath. Numerically obtained 
solutions of the Langevin Equations allow us to study the entanglement generation between the 
oscillators in terms of the time evolution of the logarithmic negativity. Our results confirm and 
extend our previously obtained findings, namely that significant entanglement between oscillators 
embedded in a free bosonic bath can only be achieved if the system are within a microscopic distance. 
We also consider the case where the bosonic spectral density is substantially modified by imposing 
boundary conditions on the bath modes. For boundary conditions corresponding to a wave-guide 
like geometry of the bath we find significantly enlarged entanglement generation. This phenomenon 
is additionally illustrated within an approximative model that allows for an analytical treatment. 

PACS numbers: 03.67.Bg, 02.50.Ga, 03.65.Yz, 03.67.Mn 



I. INTRODUCTION 

Consider two remote microscopic quantum-mechanical 
objects that are locally coupled to a common bosonic 
heat bath. It is assumed that there are no direct inter- 
actions between the microscopic objects. In this situa- 
tion any entanglement of the two objects tends to rather 
rapidly decline with time, as it is described by the general 
phenomenon of quantum-mechanical decohcrence [T], . 
On the other hand, by exchanging bath bosons the ob- 
jects do indirectly interact with each other, which, as like 
any interaction, tends to entangle the two objects. De- 
pending on which of the two effects of the bath is dom- 
inant, the two objects in an initially separable common 
state will either stay separable or develop some amount 
of entanglement. To find out whether under some given 
conditions the former or the latter case is realized is gen- 
erally a difficult theoretical problem which has been ad- 
dressed in numerous publications during the last decade 

While there seems to be agreement that for sufficiently 
small distances entanglement develops under quite gen- 
eral conditions, it is still controversial whether the mech- 
anism under discussion can also generate entanglement 
between objects in a macroscopic distance. Here we ex- 
tend and generalize our previous work on the distance de- 
pendence of entanglement generation via a bosonic heat 
bath [isj . Our main finding is that for a heat bath ex- 
hibiting a continuum of low frequency modes dissipative 
couplings will always limit the range within which entan- 
glement can be generated to very short distances. 

In Sec. In] we formulate a model system consisting of 
two harmonic oscillators that are locally and linearly cou- 
pled to a bosonic heat bath. We use a positive definite 
Lagrangian which generalizes the one that has been used 
by Unruh and Zurek in a similar context [l^ . After trans- 
forming to a Hamiltonian description we end up with a 



Hamilton operator that includes precisely the non-local 
counter term that we motivated by a different reasoning 
in 15 1 . While in [l^ we used spectral coupling density 



of states of Drude type, here we consider spectral den- 
sities J(w) with exponential cut-off at high frequencies 
and power-law behaviour J(cj) ~ uj^ with s = 1 and 3 
at low frequencies. The bath dimension is assumed to be 
either one or three. 

We derive a set of two coupled quantum-Langevin 
equation that describes the time evolution of the two 
local oscillators. From these equations we obtain the 
time evolution of the entanglement of the two oscilla- 
tors, measured in logarithmic negativity, essentially by 
use of numerical Laplace transformation. Here we cm- 
ploy an efficient new numerical method jl7j that allows 
for a better exploration of parameter space (Sec. IIIIBp . 

Strongly non-resonant frequency modes are the main 
source of decoherence, while modes with frequency near 
oscillator (spin) frequency fio tend to entangle the two 
distant oscillators. Hence, suppressing the former and at 
the same time enhancing the latter should increase the 
entanglement of the oscillators. An extreme case of this 
type is obviously realized by a single-mode cavity that is 
coupled to two systems with resonant energy levels. It 
is known that in this system entanglement can be effi- 
ciently generated. A variant of such a system that might 
be useful for entanglement over wider distances is a tube, 
i.e. a long cavity. By the boundary conditions imposed 
on transversal modes, here the spectrum exhibits a gap 
ranging from zero frequency to some finite Wg > 0. We 
analyze the entanglement generation of these type of sys- 
tems first with the numerical methods that we already 
used before in Sec. IIVI and then analytically within a 
simplified model system in Sec. |Vl In comparison to a 
free heat bath, we find for the tube significant enlarged 
entanglement. We conclude with a summary and a dis- 
cussion of our results in Sec. I VII 
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II. MODEL 
A. Lagrangian 

The model within we investigate bath mediated gener- 
ation of entanglement consists of two identical oscillators 
of frequency fJo located at positions ±r/2 and locally 
coupled to an otherwise free scalar field (^(x). The two 
oscillators represent the two remote microscopic system 
and are referred to as system oscillators in the following. 
Their canonical variables are denoted by Q1/2 and P1/2, 
their masses are set to unity. In order to avoid runaway 
solutions we choose a bilinear coupling to the velocity of 
the scalar field [l^ and thus arrive at a Lagrangian 



coupling of the system oscillators. It is merely a simple 
consequence of transforming from the Lagrangian formal- 
ism (in which the model evidently is local, cf. ffl) to the 
Hamiltonian formalism. In our previous work [15| we di- 
rectly formulated a model in the Hamiltonian formalism, 
which made it necessary to add by hand an appropriate 
counter-term in order to ensure stability and causality of 
the resulting dynamics. In fact, the resulting Hamilto- 
nian in [isj agrees with the present one up to an irrelevant 
change of oscillator variables Qi Pi. 



B. Equations of motion 

It is a straightforward task to derive the Heisenberg 
equations of motion for the environmental modes and 
the system oscillators. 



d^x0(x)(.9(x-- 



•(1) 



The coupling function (?(x) is assumed to be peaked 
around and to have a width of order 5x, describing 
the minimal distance which can be resolved by the bath 
modes 0(x). 

Quantizing the model leads to the positive definite 
Hamiltonian 
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Confining the field within a large D-dimensional box of 
size I 3> r and imposing periodic boundary conditions 
the field operators (/)(x), II(x), and the coupling function 
(/(x) can be decomposed in Fourier components as 



0(x) = 


1 

lD/2 
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n0(x) - 
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lD/2 
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.9(x) = 
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Expressed with these modes, the Hamiltonian reads 



Ql/2 = Pl/2, 

A/2 = + E-9-knke 



±*kr 



(3) 



(4) 



E-9k.9-k((3i/2 + Q2/ie''"'), 

k 

gi + e^'"-g2), (5) 

(6) 



^k = Ilk - .gk e 



^kr. 



Hk = -fc20k 



In our analysis presented below we always assume that 
at time t = the system oscillators are prepared in some 
initial state, independent from the environmental state. 
For t > both system oscillators and environment evolve 
according to the Lagrangian Eq. ([T]) or equivalently ac- 
cording to Eqs. ([3|) to (|6]). For t < the environment is 
assumed to evolve freely, without coupling to the system 
oscillators. Wc therefore require (/>k(i) and 7rk(t) to be 
solutions of Eqs. ([3]) to ([6]) for t > and to agree with 



0k, horn (0 



for t < 0. Here al^ and Sk denote creation and annihila- 
tion operators of the field mode of wavevector k. 

It can be easily checked that these field operators are 
given by 



Eff-kHk (e*''""Qi + e-i'^"-g2 

k 

4E5-k.9k (0? + 2Qig2e^'^M 



(2) 



The appearance of the operator product Q1Q2 in the 
last term of the Hamiltonian does not indicate a direct 



<Ak(t) = 0k.hom(t) 

-gkj^ dt' COs{k{t - t')) (e-5krQ^(i') + e5krQ2(i')) 

nk(t) = nk,hom(<) (7) 
+gk^ d<'fcsin(fc(t-t')) (e-^'^'"Qi(0+e*'''"Q2(0) • 
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C. Quantum Langevin Equations 



D. Covariance matrix and logarithmic negativity 



Inserting these expressions into equations ([3]) and Q, 
one obtains two Quantum Langevin Equations (QLEs) 
for the dynamics of the system operators Qi{t) for t > 0, 

Qi{t) + V^lQ^it) + dt' [ ro(t - t')Q,{t') (8) 

+ Tr{t~t')Q-,{t')]^B,[t), 

where = (1,2) and (2,1). Here, we introduced a 

damping kernel 

^(r,^) = ^gkg-kCos(fc^)e*^ 

k 

and generahzed forces given by bath operators 
Bi/2{t) = 5l3-knk,hom(0e^'^ • 

k 

It is worth emphasizing that these bath operators evolve 
freely in time; the back-action of the two system oscil- 
lators is solely contained in the memory terms of Eqs. 

®- 

The two QLEs can be written in a more convenient 
form if we combine the variables of the system oscillators 
in a vector 



m = (Qi,Q2,A,A)^, 

define a generalized force vector by 

B(t) = (0,0,Bi(t),B2(i)), 
a mass-frequency matrix 

0-10 
-1 
nl 

nl 



and, finally, a matrix 
C(t) = 



/ 





o^ 





r(o,i) r(r,i) 

\T{r,t) T{0,t) Oy 
With these definitions the QLE states 



(9) 
(10) 

(11) 
(12) 

(13) 



f{t)+Zy{t) 



d 
It 



dt'C{t - t')y{t') = B{t) . (14) 



Formally, its solution y{t) for initial y(0) and force B(t) 



m - G{t)m 



dt'g{t - t')B{t') . 



(15) 



Correlations and entanglement of the system oscilla- 
tors can be studied on the basis of the covariance matrix 

C0Vl„^it) = Trs[{yi{t),yrnit)}psm , 

where ps (t) denotes the density matrix of the system os- 
cillators. Demanding that system and environment are 
initially in a factorizing state, /5sb(0) = /5s(0) (S) /5b, the 
time evolution of the covariance matrix can be expressed 
via Eq. (|15p and a bath correlation matrix 



/C(t)=TrB[{B(t),B(0)t}/5B] 



o o^ 



KiO,t) K{r,t) 

J<{r,t) K{0,t) Oy 



as 



Cov{t) = g{t)Cov{o)g{ty (16) 
+ f dt' f dt"g{t' - t)/c(<' - t")g{t - i")^ . 

Jo JQ 

Throughout the paper we will restrict ourselves to the 
consideration of Gaussian states which are entirely de- 
termined by their covariance matrix. Then, a convenient 
measure for the entanglement is the logarithmic negativ- 
ity 

Sa/- = - ^ log2 (min(l, Ai)) , 

1=1,2 

where X1/2 are the symplcctic eigenvalues of the partial 
time-reversed covariance matrix Cov'^^ = PCovP with 
P = diag(l, 1, 1,-1). Loosely speaking, the logarithmic 
negativity measures the deviation of the oscillators being 
in a separable state p^ . 



E. Spectral coupling density of states, thermal 
bath correlations 



The model needs to be further specified by fixing the 
oscillator-bath couplings gk- To this end we assume 
isotropy, i.e. = gk, and define a spectral coupling 
density as 



J(c^) =^|.gkpc^5(w-fc) 



which in turn we suppose to be of the form 



(17) 



(18) 



where g{t) is Green's function of the QLE p^ . 



Here, 7 denotes the overall coupling strength, s is the 
spectral index, and flc is a cutoff frequency. Note that the 
definition of the spectral coupling density differs from the 
usual definition in the extra oj-factor, which attributes to 



4 



the fact that we consider a coupling to the field velocities 
rather than to the field itself. 

With this definitions we obtain for the damping kernels 
the explicit expressions 



r(r,i) 



dui ^ cos(a;t) cos(ajr) 

'd^^!fMcosM^^^ 



(19) 



for a one-dimension and a three-dimensional bath, re- 
spectively. When the environment initially is in a ther- 
mal state pb of temperature T, the bath correlations 
determines to 



K{r,t) 



dujJ^^iuj) coth (^) cos(a;t) cos(a;r) 
dcoJ^^ica) coth (#) cosM)^ . ^ ' 



III. NUMERICAL SOLUTIONS 




FIG. 1. Logarithmic negativity E as function of time (in units 
of l/f2o) for different distances for the ID bath (in gray) and 
the 3D bath (in black). Dashed lines represent asymptotic 
values. The parameters are the same as in Fig. O A zoomed 
in view of the beginning of the evolution t < 0.6/f2o is shown 
in the inset. 



Methods 



The solution of the integro-differential equations ([5]) 
which is formally given by the expression was ob- 
tained numerically. 

In order to solve the homogeneous part (jS]) in the po- 
sition basis we introduce the auxiliary functions 



Mr,t,u)= / dt' {r{0, t - t' + u)Q,{t') 
Jo 

+Tiv,t-t' + u)Qj{t')) 
which satisfy the partial differential equations 
dtMr, t, u) = r(0, u)Q,{t) + r(r, u)Qj{t) + a„/,(r, t, u) 
By means of the Fourier transforms 
1 



Mr,t,u) = 



r(r,u) 



/27r J- 



1 



d.se'^'Mr,t,s) 



ds e*""r(r,s) 



we obtain the ordinary differential equations 

9t/,(r, t, s) - f(0, s)Q,{t) + f(r, ,s)Qj(t) 

+isf,ir,t,s). (21) 

Discretizing the inverse Fourier transform of f{t, s) we 
can rewrite the homogeneous part of ([8]) according to 



2rig,-id + l 



k=l 



dth{v,t,sk). (22) 



where As^ denotes the grid spacing. Equations (PT|) and 
((22)) form a coupled system of ordinary differential equa- 
tions which can be solved by standard algorithms. 

Using the numerical solutions for and Pj . the Greens 
function can be deduced from the homogeneous part of 



(1151) . Futhermore, in order to obtain the covariance ma- 
trix ([T5| . three numerical integrations have to be nested, 
two for the time integrations and one for the bath cor- 
relator (PH]) which is not available in closed analytical 
form. 



B. Free Space Environment 

We now present our results from the exact numeri- 
cal solution of our model. For the free bath, which is 
discussed in the following section, we used a grid size 
of rigrid = 10000 and constant grid spacing Asfc with 
Smax = (2ngrid + l)Asfe = lOfic. As wc want to study 
the creation of entanglement between the two system os- 
cillators, we start with a separable state ps(0). For sim- 
plicity, ps{Q) will always be the ground state in the fol- 
lowing. The bath is initialized in a thermal state pT with 
zero temperature. 

We look at both the asymptotic values reached af- 
ter the system has equilibrated, as well as the dynamics 
for short and intermediate times. The free parameters, 
which remain after scaling out mass and frequency of 
the oscillators, are distance r, cut-off frequency Oc, and 
damping 7. Generally, we measure distances in units of 
c/r^o and frequencies in units of 

The full time evolution exhibits fast dynamics at the 
beginning due to the coupling of the bath being switched 
on suddenly, see Fig.[T] Later the entanglement oscillates 
until it reaches its asymptotic value. For the physically 
unrealistic case r = 0, the relative position Q2 — Qi re- 
mains undamped during the whole time evolution. The 
damping increases for larger distances r. From Fig. [1] 
it can be deduced that the asymptotic value is reached 
more quickly for larger distances. 

As the maximum entanglement reached at short times 
is larger than the asymptotic value, one can ask the ques- 
tion if it features a more favorable distance dependence. 
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FIG. 2. Smallest distance rmax at which the oscillators stay 
separable during the entire time evolution as a function of 
the inverse cut-off 1/ilc for the ID bath (in gray) and the 3D 
bath (in black). The parameters are the same as in Fig. [31 
The inset shows rmax as function of the coupling strength 7. 



Wc have plotted the smallest distance at which the oscil- 
lators stay separable during the entire time evolution as 
a function of the inverse cut-off in Fig. [5] While the dis- 
tance is somewhat larger than in the asymptotic case, see 
Fig. m it can still be upper bounded by a linear function 
proportional to l/flc- Therefore the qualitative behav- 
ior is the same. Exactly as in the asymptotic case, the 
coupling strength 7 only has a minor influence on rmax- 
The asymptotic state can be calculated by taking the 
Fourier transform of the Langevin equation (|14p [l9| . 
from which we have dropped the quickly decaying term 
y(0)C{t). The resulting Quantum Langevin equation 

y{t) + Zy{t) + f dt' e{t ~ t')C{t - t')y{t') = B(t) 



is transformed into the algebraic equation 

— icjy(w) + Zy{ijj) — iujTZ{uj)y{uj) = B(w) 
with the matrix TZ{lu) defined by 

C{uj) = - / dte"^'C{t) 

^ J -QO 

Tliuj) = / dt e'^'C{t) = C{u) + -P / dw'-^ . 
Jo ^ J -00 w'-w 

The Fourier transformation of C{t) cancels with the in- 
tegration over uj, so that the matrix elements T{v,ljj) of 
C{uj) corresponding to the elements r(r, t) of C{t) are 
simply given by 



f(r,^) 



4|w| 



cos(a;r) . 



The equilibrium equal time correlation function of the 
system can then be expressed in terms of the matrix 



0.06 



0.04 - 



0.02 - 




0.005 0.01 0.015 0.02 0.025 0.03 

r[c/r!o] 



FIG. 3. Asymptotic entanglement of the system oscillators 
measured in logarithmic negativity i5 as a function of distance 
r (in units of c/cjo) for the ID bath (in gray) and the 3D bath 
(in black), damping constant 7 = f^o, and cut-off frequency 
Sic = 10^0- E drops to zero at a rather small critical distance 

' max 

< c/fic (cf. Fig.©. 
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l/0c[l/0o] 



0.1 



FIG. 4. Distance rmax at which the oscillators become asymp- 
totically separable for the ID bath (in gray) and the 3D bath 
(in black). The parameters are the same as in Fig. [S] 



as the integral 



({y.(0),y,(0)}) = — • 

"0 



E 

fc,;G{3,4} 



du:Fik{Lo)Fji{-bj)J{\uj\)cos{u:r5ki) ■ 



The behavior of the asymptotic entanglement as a func- 
tion of distance r is plotted in Fig. [3] It drops to zero for 
rather small distances rmax, which arc on the order of the 
inverse cut-off l/f2c as shown in Fig.|4l This is only very 
weakly dependent on the coupling strength, which can 
be explained by the fact that 7 controls both coupling 
and decoherence at the same time. Changing it cannot 
be used to increase the entanglement capabilities of the 
bath. 
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IV. WAVEGUIDE ENVIRONMENT 

While environmentally created entanglement does not 
seem to be possible in free space, changing the geome- 
try can have a big impact. A simple idea is to make the 
bath effectively one-dimensional by placing the two os- 
cillators inside a waveguide with quadratic cross-section 
of side length a. The quantization in transverse direc- 
tion suppresses all frequencies below the first mode with 
frequency wq = V^tt /a. In order to keep things simple, 
we will for now only consider this first mode in trans- 
verse direction. The spectrum in longitudinal direction 
will remain continuous. 

We split the wave vector k into its components 

27rrii 2tt 

k = — — ei H (71362 + 71363) . 

L a 

with ni e Z \ {0}. 

The component in ei-dircction can approximated by 
an integral due to the large length of the tube in this 
direction. For this, we make a change of variables ni — > 
|k| = fc which implies 



dni 



I 



k dk 



27r 



A:2 



At the same time we will examine what happens to an ad- 
ditional term e**^"" as it occurs in the damping kernel (fT9|) 
and in the bath correlator ([20| . Assuming that r || ei we 
find for the damping kernel 



r(r,t)= / dc^V|5k|'<5(w-fc)e^'''" 
Jo 



(23) 




0.5 LOo 



FIG. 5. Coupling spectral density J (a;) of a quasi one- 
dimensional bath. We only consider the first transversal ex- 
citation and approximate the rest of the spectrum with the 
free bath (black line). The full spectrum is shown in gray. 



with 

Jwg('i^) 



87a;, 



(24) 



Similarly, the bath correlator reads 
i^wg(r,i) = 

j duo Jwg(w) coth cos(cj<) cos{r^ — ujq) . 

In order to approximate the full spectrum, we will use 
the approach depicted in Fig. [SJ Only the first excitation 
is considered and the rest of the spectrum is replaced 
with the free bath. This is a good approximation as long 
as the oscillator frequency fio is close to the frequency of 
the first excitation wq. 



du 



E 

"2,113=1 



cos ( rv^w2 - ^(n2 +7i§; 



where the factor 8 is the number of octants of a sphere. 
Taking in the limit — !• and comparing with 
we find 

2 _ J{^)^1 

where we assumed isotropy of the spectral couplings, that 
is gk = 9k- If the oscillator frequency fig is close to the 
first excitation, we can approximate the sum in (|23p by 
the first term 712 = 713 = 1. Thus, the damping kernel in 
the waveguide reads 



Fwg(r,i) 



' dw^ 



cos(wi) cos (t^J' 



V. NUMERICS 

Due to the inverse square root divergence of the cou- 
pling spectral density, the required numerical effort is 
larger than in the free case. We were able to achieve 
good convergence by using the non-uniform distribution 
of grid points. 



-^max ^0 7 2 , 

k + UJQ 



("grid) 2 

resulting in the fc-dependent spacing 
2k 



Ask = 



(%rid)' 



The results of the numerical calculation are shown in 
the form of a density plot in Fig. |6l Shades of gray encode 
the value of the entanglement which is drawn as function 
of time (x-axis) and separation (y-axis). 

We again start with a separable state, which now is a 
single mode squeezed state with the squeezing parameter 
K = 10 

Cov(O) = diag(l/(Kr2o), l/C^rJo), nflo, Kflo) . (25) 
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FIG. 6. Dependence of the time evolution of the entanglement 
Ej\f on the spatial separation r if the systems is placed inside 
of a waveguide. The value of the entanglement is encoded in 
the shade of gray as indicated by the scale on the right side. 
Additionally contour lines are drawn at the values Ej^ — 0, 
0.5, 1, and 1.5. The parameters are fie = lOno, 7 = 0.05 fio, 
and T = fio with an initially separable single mode squeezed 
state K = 10. 
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FIG. 7. Distance dependence of Ej^ &tt = GO/fio for different 
temperatures. All other parameters are the same as in Fig. [S] 



VI. ANALYTICAL MODEL 

As we have seen in section (|V)l . boundary conditions 
which are imposed on the bath modes have an impact on 
the generation of entanglement. According to equation 
([24)) . the spectral weight of the frequencies in an interval 
5u) centered around the van Hove peak is significantly 
larger than the spectral weight of frequencies in an in- 
terval 5u which is distant from the singularity. This fact 
can be used for a derivation of approximate analytical 
solutions of the model discussed in the previous section. 
Therefore, we take into account only the first van Hove 
singularity at = ^/Sn/a and replace the coupling spec- 
tral density ((24)) by two effective macroscopic oscillators, 
representing the coherently oscillating symmetric and an- 
tisymmetric field modes around the van Hove singularity. 
Furthermore, we take into account for the incoherently 
oscillating modes of the bath by a background spectral 
density, given by (jl8p . This replacement, schematically 
depicted in fig. [8l will be a suitable approximation as 
long as the modes within the van Hove peak are oscil- 
lating coherently, that is, for t,r < 1/Suj. The effective 
oscillator representing the symmetric modes will deter- 
mine strongly the generation of entanglement whereas 
dissipation and decoherence is determined by J^^ . 




Jwg{LJ,uJo) — > J^^{(jj) -h oscillators 



The ground state corresponds to k = 1. Wc note that 
increasing the squeezing increases the absolute value of 
the entanglement but does not change the distance over 
which it is generated. The other parameters are com- 
parable to the free case: The coupling strength is set to 
7 = 0.05 VIq and the cut-off frequency is fie = 10 Oq. This 
time, however, entanglement is generated over distances 
which are more than two orders of magnitude larger than 
in the free bath. The effect is most pronounced if the 
frequency of the peak loq is close to the frequency of the 
oscillators VIq. 

Fig. [7] shows the exact distance dependence. The plot 
corresponds to a vertical cut through Fig. |6] at time 
t = 60/ilo. Increasing the temperatures decreases the 
entanglement generation, but the effect persists even for 
comparably large temperatures of order fio- 



FIG. 8. Schematical illustration of the replacement of the van 
Hove coupling spectral density by the free coupling spectral 
density and effective oscillators with frequency (jjq. 

We choose r = rei and decompose the system vari- 
ables and the field modes in symmetric and antisymmet- 
ric compositions with respect to the ei-direction, that 
is 

nk,s/A ^ (nfci,fc2,fc3 ^^-kiMM^ 
Qs/A = ^iQi±Q2) 

We find that the effective coupling oscillator interacts 
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only with the symmetric mode Qs- Within this approx- 
imation, the Hamihonian ^ reads 



H — Hb^s + Hb^a + Ho,s + Hq^ 



Hi, 



with 



B,S/A 

Ho,s 

Hint 



\ E (ft 



k,5n_k,S/A 



'^\i.SV--k,S/A 



ki>0,k2.k3 



s 



2 -2 ^ 



1 / "2 , 2 »2 \ 

' 2 \Pa+^o<1a) 



gpsQs : 



All >0,fc2 ,A;3 



IIicsOs cos 



kir 



-iIIk.AQAsin ( ^ 



Here we denoted the canonical position and momentum 
variables of the effective oscillators with qs/a and Ps/Ai 
respectively. The coupling constant g depends on the 
coupling spectral density and will be determined below. 
Furthermore, the counterterms of the Hamiltonian lead 
to a renormalization of the system oscillator modes, that 
is 



cLu 



1 ± 



sin(wr) 



The Hamiltonian Hg can be diagonalized by means of 
the transformations 



Qi = 

Pi = 



Ps + £,uJoq 



P2 = 



2guJo 



(26) 



and adopts the canonical form 

1=1,2 

+ \{Pl^-n\Q\+p\+u:litA)). 
The eigenmodes are given by 



The time-evolution of the system oscillators and the cou- 
pling oscillator will be treated exactly, for the incoher- 
ently oscillating modes we will use a master equation 
approach. 



The time evolution of the density matrix is determined 
by the differential equation pss = —i[H,PsB\- Assuming 
that the bath correlators decay on a time scale which is 
much shorter than the time scale of the system oscilla- 
tors allows for a Born-Markov-approximation 0. The 
density matrix of the system factorizes in a symmetric 
and an antisymmetric part, ps = ps ® PA, which obey 
the differential equations 



dtps I A 



-i E \_^s/A,PS 

i=l,2 



(27) 



dt vs/A{t, r)[Qs/A, [Qs/A{-t),Ps/A] 



+i / dt Us/Ait, r)[Qs/A,{Qs/A{-t), ps/a}] 
Jo 

The bath correlators vs/A and Ps/A are given by 
'^s/A{t,r) 



did J ' 



3D 



, ^ , / uj \ f sin(ci;r) \ 

{lo) coth (— j [l ± -^^j cos(wi) , 



Ps/A{t, r)^\j duJ^^ico) (^1 ± Met) . 

Decoherence and dissipation is completely deter- 
mined by the double commutators involving eight 
different system-bath-corrclators for the symmetric 
and antisymmetric bath modes, respectively. The 
decoherence rate is determined by system-bath- 
correlators ai^s/Ai^) = dti's/ACOs{flt), the 
correlator Q!2,s/a(^) ~ — i/g/yi sin(r2t)/r2 is 
called anomalous-diffusion coefficients, a^g/ji{Q) = 
~ Io° dt fig/j^^cos{nt) introduces a the lamb shift and 
finally 0:4^5/^1(0) = dt ^g/j^sm{Qt)/Q determines 
the strength of dissipation. The frequency O is equal 
to f2i or ^2 for the correlators with index S and equal 
to uj or flA for the correlators with index A. Explicit 
expressions are given in the Appendix, see equations 

dMli-gl]). 



A. Induced oscillator coupling 

The effective coupling constant g can be roughly esti- 
mated from the expressions of the correlators vs or /zg. 
Splitting the wave vector as in section pvp and replacing 
the rii-summation by an integral we find 



poo 2. z • --A) 

'^s{t,r) = duj y2 

JO 1 

n2,n3 = l 



"2, "3=1 7TLij\/ uj'^ — ^{n2 + n,'^) 



X coth cos(cji) cos^ 



Restricting ourselves to the first van Hove singularity, 
71,2 = 'T-2 = 1, and integrating from loq to ujq + Scj, we 
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deduce the relation 



g^{p{t)ppB) 



9^ Y cos (wot) coth(^2^^ 



i^^j3^(c.o)J— cos(wot)coth(^) 

TT V 1^0 \ li / 



The modes within the van Hove peak are osciUating co- 
herently for distances r < cjbuj which introduces an effec- 
tive distance dependence of g. Since this coupling con- 
stant would increase indefinitely for r — > 0, we need a 
natural cutoff for the smallest distance possible which is 
given by l/r^c- Therefore we have (5a; = 1/(7' + l/f^c) and 



a = .9(0 



2J3^(wo) 



Wo?' - 



(28) 



B. Analytical solutions 

Although we restrict ourselves to Gaussian density ma- 
trices, equations (P7| lead to a system of coupled nonlin- 
ear differential equations in the position representation. 
For this reason we transform the density matrix into the 
"fc — A" -representation [l^. Since the procedure is the 
same for ps and p^, we restrict ourselves to the former 
one. The transformed density matrix is given by 



kQ+AP 



p(k, A) Tr pse 



with the vectors k = {kxMY , ^ = (Ai,A2)^, Q = 
(Oi; O2) and P = {Px.P-i). This representation is related 
to the Wigner-Distribution via a double Fourier transfor- 
mation and has in position representation the form 



p(k, A) = / dxxdxie^^^ps x 



A 



A 

T 



where x = (x\,X2) labels the diagonal elements of the 
density matrix. Using this particular representation we 
find that the master equation is linear in the derivatives 
with respect to hi and A^. With an Gaussian ansatz for 
the density matrix /^(k, A) we end up with a linear first 
order system of differential equations, see equation ([34|) 
in the Appendix. 

The dissipative part of equation (P7)) leads to a cou- 
pling among 14 differential equations, thus the eigen- 
modes cannot be found analytically. However, neglecting 
terms 0(7^) in the dissipative part of equation p9|) al- 
lows for an analytical treatment of (j34p . After removing 
the higher order couplings we are left with the usual 4 
double commutators known from the Caldeira-Leggett 
model (see e.g.Q). 

The expectation values of the anticommutators, and 
therefore the negativity, can be expressed in terms of 
the solutions Ci (t) of the differential equations (see equa- 
tions [35]) . From the approximate solutions ((36|) - p7p we 
find that the damping of the eigenmodes is completely 



determined by the dissipation-correlator. The lambshift 
modifies the eigenmodes fti/2- We found that analytical 
expressions coincide very well with the numerical inte- 
gration of the differential equations p4p for small and 
intermediate times. 
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FIG. 9. Dependence of the time evolution of the entanglement 
on different squeezing parameters. The remaining parameters 
are s = 3, 7 = O.05f2o, = 3f2o and r = 2[c/Qo]- 
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FIG. 10. Dependence of the time evolution of the entangle- 
ment on different temperatures. The squeezing parameter is 
K = 5 and the remaining parameters are the same as in Fig. [9] 



C. Results 

We consider the system oscillators to be initially in a 
Gaussian state with squeezing parameter k, see equation 
ipS]) . the effective coupling oscillator is assumed to be 
in the ground state. From Fig. IH] we deduce that the 
maximum value of the negativity depends crucially on 
the initial squeezing. Due to the dissipative time evo- 
lution, the system oscillators relax within intermediate 
times from the squeezed and nonentangled state to an 
entangled state with lower energy. In contrast, preparing 
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t[l/Oo] 

FIG. 11. Dependence of the time evolution of the entangle- 
ment on the spatial separation of the system oscillators with- 
out coupling oscillator. The parameters are k = 5, s = 3, 
7 = 0.05^0, fie = 3^0 and T = 0.01 [^o]. 
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FIG. 12. Dependence of the time evolution of the entangle- 
ment on the spatial separation of the system oscillators. The 
parameters are k = 5, s = 3, 7 = 0.05Slo, ~ 3fio and 
r = 0.01[f7o]. 



the system initially in the ground state, that is k — 1, 
only negativity of order 7^ /rig generated. 

The temperature dependence is very weak for T ^ rig 
whereas for T = O(r2o), the negativity decreases signifi- 
cantly, see Fig. [TOl 

In virtue of the van Hove singularity we find significant 
entanglement over distances r 0{c/ilc) as can be de- 
duced from Fig.[T2l whereas for g = 0, which corresponds 
to the time evolution of the free space environment, we 
find no entanglement for distances that are larger than 
0{c/nc), see Fig.HH From the expressions ^ and ((55)) 
we see that the generation of entanglement for 17 = is 
related to the difference between the dccohcrence and dis- 
sipation rates of the symmetric and antisymmetric mode, 
respectively. 

Although various approximations were necessary in or- 
der to derive analytical expressions for the time evolu- 
tion, we find good qualitative agreements with the exact 
numerical calculations, compare Fig. [T^] with Fig. [SI 

Taken the approximate relation (pS)) for granted, all 



entanglement vanishes in the asymptotic limit. How- 
ever, treating g and 7 independently from each other 
we find for sufficiently large g nonzero asymptotic neg- 
ativity which can also be determined from the thermal 
expectation value of the covariance matrix up to correc- 
tions of order 7^. In Fig. 1131 we depict the temperature- 
dependence of the asymptotic entanglement for g = 
for different values of r. From Fig. [T3] we deduce that a 
critical minimal value of g is necessary in order to observe 
asymptotically a nonzero negativity. 
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FIG. 13. Dependence of the asymptotic entanglement on the 
temperature for different spatial separation of the system os- 
cillators. The parameters are s = 1, 7 = O.Olfio, fic = SQo 
and g = l[f2o]- 
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FIG. 14. Dependence of the asymptotic entanglement on the 
coupling g for different temperatures. The parameters are 
s = 1, 7 = O.Olf^o, f^c = 3f7o and r = l[c/r2o]. 



VII. CONCLUSIONS 

We investigated the generation of entanglement of re- 
mote quantum systems via a bosonic heat bath. Starting 
from an Lagrangian formulation describing the coupling 
of two remote harmonic oscillators to a scalar field envi- 
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ronmcnt, wc derived the Langevin equations of motion. 
The coupHng of the oscillators to the field momenta guar- 
anteed the positivity of the corresponding Hamiltonian. 
This led to an counter-term which was already proposed 
in [l^ in order to preserve causality and to avoid runaway 
solutions. 

In case of a free space environment the generated en- 
tanglement drops quickly to zero if the distance between 
the ocsillators exceeds the spatial extension of the indi- 
vidual quantum systems. For one and three spatial di- 
mensions, we found that this statement is independent 
of the initial conditions. 

In contrary, by imposing boundary conditions on the 
heat bath it is possible to generate entanglement over dis- 
tances which exceeds significantly the spatial extension 
of the quantum systems. In particular, we considered 
a waveguide with a transversal extension corresponding 
to the inverse system frequency. The boundary condi- 
tion leads to a van Hove peak in the coupling spectral 
density. The entanglement can be enhanced significantly 
if the van Hove peak is in resonance with the remote 
quantum systems. Furthermore, the situation can be im- 
proved if the quantum systems are prepared initially in 
a strongly squeezed state. 

Considering the coherently oscillating modes as an ef- 
fective oscillator coupling the quantum systems to each 
other, we derived analytical solutions which are qualita- 
tively in good agreement to the numerical results. 



erators (3i/2 and P1/2, that is 



^11 



Pf + n. 



2r\2 



1=1,2 

ai(f^i) 



(29) 



^_^^2 > lQiJQi,P]] + ^[^2,[Qi,P]] 



^2 



-^^^ (ah, [Q2,p]] + n^iQi, [Q2,p]] 



[Oi,{Oi,p}] + ^[A,{Qi,p}] 



i+e V'"" ' ^2 



[Oi,{A,/5}] + i-[P2,{A,/5}] 



i + e v^------ ■ 



ii+e)^'2 



aP2AP2,p}]+^2[Ql,{P2,p}] 



Xai{^2) 



1+e 



aP2,{Q2,p}]+^2[Qu{Q2,p}'_ 



VIII. APPENDIX 



where f was defined in equation (j26|) and the a, are corre- 
lation functions that are given explicitly in section lVIII CI 
The equations of motion in the "/c — A" -representation 
can be constructed by wrapping 



A. Differential equations 

In order to avoid a cluttering of indices we will restrict 
ourselves to the density matrix for the symmetric mode 
and suppress the label 5". The equations determining the 
antisymmetric mode follow from a replacement of the 
eigenmodes and canonical variables in the limit ^ — > 0. 
The time evolution for the oscillator variables is given by 



Tr LD... = tr 



I kQ+AP 



over equation ([32]) . Using the relations 



(30) 



Ql/2it) = Qi/2COs(rJi/2t) 



ill 



/2 



-Pi/2sin(f7i/2t) , 



whereas for momentum variables we find 



Pi/2it) = Pi/2 cos(rJi/2i) - fli/2(5l/2Sin(rJi/2t) 



Using relations (1^51) we can decompose the differential 
equation (|27p with respect to the (time-independent) op- 



Tr(i?[g„[Q„p]]) = A,Ajp 
Ti:{D[Q,,[P,,p]]) = -A,kjp 

Tr(^[A,[4-,/5]])--^A/5 

Tr(D[4,{4,/5}]) = -2iAA^.p 
Tr(£)[4,{A,P}]) = -2iA,aA,p 
Tr(^[A,{4,P}]) = 2zfcA,p 
TTiD[h{P,,p}])^2ihdA,P, 



(31) 
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we arrive at 



2$ 



Y^«3(f^l) ( Ai - 

2^ 



dkiP 



1+e 



1+e 



2a4(r!2)(f)2Ai-^fc2) ) dkj 
-ai(ni) fA?-^fc2Ai 



^2 



-a2{^i) ( fciAi - 



ai(n2)(f^2A:2Ai-^A:|) 



+$a2(f^2)(f^2AiA2-Cfc2A2) 



p. (32) 



With the Gaussian ansatz 



PQs,q = exp [ - cikl - cafciAi - c^Aj - ic^ki - icgAi 

-Cgfcj - C7fc2A2 - C8A2 - iCgk2 - 10x0^2 

-ciifcifca - ciafciAa - ci3fc2Ai - ci4AiA2](33) 



C6 



;a4{n2)c6 + 1 + 



2e 



1 /2e 



(l+^2)f^2 

2e 



a3(f^2) C7 



1 + V^^: 

+ ^ai(1^2) 
-2 2^2 

C7 = -2f^2C6 - 

4^2 



^a3(fJi)cii + --a4(rJi)ci3 



^2 



a4(Jl2)c7 



2 + 



2^ 



2^ 



^a3(fJi)ci2 + ^a4(rii)ci4 + e a2{^2) 



cg = Clo + 1—7-72 ( ^"3(i^l)c4 + ^Q!4(fil)c5 



lis 



2^2 



-2^ a4(fi2)c9 + -^Q!3(f^2)ci0 
"2 



1 / 4^ - 2^ - 

cii = , , ( :^a3(rii)ci + =-Q;4(fii)c2 
h "2 



V^^2 



2^2 



~Ci3 



04(172)211 
f^2(l+e^ 



1 



2^2 



C12 



-02(^1) 



C12 = — ri'cii + ci4 



we arrive at the first order system 



Cl = C2 

C2 = - ( 21]^ 



C3 



l + e' 



Q;3(rii) Cl + 2C3 



(^2a4(f^i)c2 - 2^1^204(1^2)011 
2^ - 

'o~'^3(f^2)ci2 + a2{^l] 
"2 



1 



i + e 



;a3{Cli) C2 



1+e 

2C 



4a4(17i)c3 - 2^n2a4{n2)ci3 



0:3(^^2)214 - ai{Qi] 



C4 = C5 



C5 



ra3(rii) C4 



^ ( 2q;4(11i)c5 - 2^1^204(02)09 



i+e\ 

-^03(02)010 

"2 



1 



Cl3 



2^ - 4^ - - - 

2 > ?r'^3(Oi)c2 + — 04(r2i)c3 + 4^17204(02)25 



2e 



ii+e)^2 

2 



03(02)27 



1 + C 



:a3(Oi) + 02 cii 



1 + e 



O4(0i) +e^a4(02)jci3 
-%To2"3(^2))ci4 



(i+eM 



(i+e')02 
1 



Ci4 



ai(Oi) + 01(02) 
4C 



2^0204(02)27 - —03(02)28 
2 



1 + e' V ^ — 17 



0^ 



1 + e' 



03(01) C12 - O2C13 



1+e 



(204(01)014+^^^202(02)) 



(34) 
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The expectation values of the anticommutators can be 
given in terms of the functions q (t) according to 



({Qi,Qi}> 


= 2(2ci 


+ cD 


({4,4}) 


= 2(cii - 


f C4C9) 


({4,4}) 


= 2(2c6 




({4, A}) 


= 2(c2 4 


C4C5) 


({4, A}) 


= 2(c7 4 


CgCio) 


({4, A}) 


= 2(C12 - 


+- C4C10) 


({A, A}) 


= 2(2c3 




({A, A}) 


= 2(2c8 




({A, A}) 


= 2(ci4- 


+- C5C10) 


({4, A}) 


= 2(ci3 - 


f C5C9) . 



(35) 



B. Approximate Solutions of the Differential 
Equations 



From the initial condition (|25p we find that the func- 
tions C4,C5,cg and cio are vanishing since they describe 
momentum and position displacements that are absent 
in symmetric Gaussians. At t = 0, the non-vanishing 
anticommutators have the expectation values 



({Ps/A,P5/a}) = 

{{Psia.Psia}) 



1 



{{^siaAsiA) 
1 

({Qs/A, Qsia)) 



= Wo 



Kilo- 



Neglecting terms of order 9^7 in the differential equation, 
some of the equations ([M]) decouple from each other. For 
the coefficients concerning the oscillator with the variable 
Qi, we find 



Cl 



3 



_A,t , ai(r^i) - 2a2(^i)Q^4(^i) 
4«4(l^i)(2a3(ni) + 172) 



C2 =E^«^*e^'* 
i=l 

C3 = E f T + ^1 + 2a3(f^i) + A,e^' 



«i(»i) 
4a4(ni) ' 

with the eigenmodes 

Ai = -2a4(17i) 



(36) 



A2,3 = -2 ± iJ^i + 2a3(f^i) - ali^i) 



The are chosen such that 



ci(0) 

C2(0) 
C3(0) 



4Kflo'^o(l + ^2) 

C2(0) =0 

Kfio + ^^Wo 



4(1 + a ■ 
For the oscillator Q2 we find 

eg = Bie^^^^t 4. ^2e-2'^"^2t 4. ^3 

C7 = 2ifl2^ie^'^"* - 2ifl2-B2e-2*^"^^* 
Cg = -fl^(Bie2'^^* + ^26-2*^^=^* - B3) . 

The -Bi are chosen such that 

Wo -I- ^^Atfio 



C6(0) = 
C7(0) 
C8(0) 



4172(1+^2) 

C7(0) = 

n^(K»o +C^Wo) 
4Kr2oWo(l + ^^) 



The coupling between Qi and Q2 is described by C11...C14, 
which read 



Cll 



C12 



2f2|(a4(rJi) + Ki) 



(37) 



4 

ci3 = E 

217|(a4(l^l) + K,;) 



+E^> 



Ci4 



E^^ 



n2(2a3(l^i)+n2-172-K2; 



2a3{ni) + ni-ni + 2a4{ni)K, + 

with the eigenmodes 

«^i,2,3,4 = —a4{Cli) ± i ^2Q;3(r2i) — a^i^i) + Clf + fl^ 

1/2 



± 2f]2V2a3(^^i) - aK^^i) + f^? 



The Ci are chosen such that 



cii(O) = 



Cl2(0) 
C13(0) 



C^2(W0 - hQq) 

2Kf}oWo(l+e2) 

^(WO ~ Kflg) 



Cl4(0) =0. 
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C. Bath Correlators 



with 



For the decoherence rates we find 
a,^s,m = \j^-m (l±^)coth(|)(38) 
The dissipation correlators read 



«4,5M(f2) = ^^(l±^^) . (39) 



The anomalous diffusion correlators have the form 



- 27r(i + .) fn 



cos(7rs)5R <^ r 



sin(7r,s)3 < T 



n 
a: 



coth I ^ I X 



8T-/ncT{l + s) 1 f ^ 



s-2 



± - f V coth f 1 X 



2 cos(rJr)e-^ sin(7rs)5R <^ T -s + 2, 



+ cos(7rs)Q <^ r -s + 2, 



a: 



K^--)r(-. + 2,o(i^-.. 



-.s + 2,-0 I ^ + 



47rr(s- 1) ^ /27rnTY^^ d„ 



Cn = e 2 1 I — s, —I- ' 



17e J 

27rnT\ 



+ e-'2-nT(J--^r)-^¥p _^ 2, -i2™r - i,' 



The Lamb shifts due to the bath correlators are given by 



a3,s/A(f^) 



27r(3 + s)i7 ( n 



e'^ i cos(7rs)3fi <| T ( -2 - s, 



sin(7rs)S F ( — 2 — s, — — 



+e^Ti-2-s,- 



47l73F(3 + s) 
7:n^{2 + s) 



7F(i + s) fn_ 

fir 



s-l 



2 cos(rJr)e~^ sin(7rs)5R <^ F -s, - 



1^ 



- cos(7rs)3 F ( — s, — ^ 



+,e"(^+'^^)F('-s,n('^+zr 



_ie"(Tfc-"-)r ( -5,17 ( — - ir 



1 



(41) 



The Lamb shifts originating from the counter terms are 
determined by 



duj 



UJ V 



± 



sin(ci;r) \ SjQc 



Q u; \ (jjr J n 



F(s) (42) 



(40) 



± 



(1 + nlr^)- — T{s - 1) sin [{s - 1) arctan(»cr)] ' 
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